Draft version February 5, 2008 

Preprint typeset using 1^1^^ style cmulatcapj v. ^jllj^A, 



in 

O 

o 

< 



> 
OO 

\o 

(N 
00 

o 
o 

6 



CMB COMPONENT SEPARATION BY PARAMETER ESTIMATION 

H. K. Eriksen\ C. Dickinson^, C. R. Lawrence'*, C. Baccigalupi'', A. J. Banday"', K. M. Gorski", F. K. Hansen'^, 
P. B. LiLJE*^, E. PiERPAOLr', M. D. Seiffert^", K. M. Smith" and K. Vanderlinde^^ 

Draft version February 5, 2008 

ABSTRACT 

We propose a solution to the CMB component separation problem based on standard parameter 
estimation techniques. We assume a parametric spectral model for each signal component, and fit the 
corresponding parameters pixel by pixel in a two-stage process. First we fit for the full parameter 
set (e.g., component amplitudes and spectral indices) in low-resolution and high signal-to-noise ratio 
maps using MCMC, obtaining both best-fit values for each parameter, and the associated uncertainty. 
The goodness-of-fit is evaluated by a statistic. Then we fix all non-linear parameters at their low- 
resolution best-fit values, and solve analytically for high-resolution component amplitude maps. This 
likelihood approach has many advantages: The fitted model may be chosen freely, and the method is 
therefore completely general; all assumptions are transparent; no restrictions on spatial variations of 
foreground properties are imposed; the results may be rigorously monitored by goodness-of-fit tests; 
and, most importantly, we obtain reliable error estimates on all estimated quantities. We apply the 
method to simulated Planck and six-year WMAP data based on realistic models, and show that 
separation at the /xK level is indeed possible in these cases. We also outline how the foreground 
uncertainties may be rigorously propagated through to the CMB power spectrum and cosmological 
parameters using a Gibbs sampling technique. 

Subject headings: cosmic microwave background — cosmology: observations — methods: numerical 



1. INTRODUCTION 

As experimental techniques improve rapidly, and new 
high-sensitivity ground-based, balloon-borne, and satel- 
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lite missions are being planned and implemented, the 
main problem in CMB measurement has changed from 
instrumental noise to separation of the cosmological 
CMB signal from non-cosmological foreground signals. 
This problem will become even more important as our fo- 
cus shifts from observations of temperature anisotropics 
to polarization anisotropies; while a simple template- 
fitting approach proved adequ ate for the first-year 
WMA P analysis i Bennett et"aLll2003albl: iHinshaw etUI 
12003(1 ■ no such hopes can be held for future polarization 
experiments. 

While component separation is a difficult problem, 
it is not intractable. Since the cosmological CMB 
radiation follows a virtually perfect black body spec- 
trum (Mather et al. 1999), whereas all known non- 
cosmological signals have non-thermal spectra, it should 
be possible to disentangle the various contributions us- 
ing spectral information. This fact motivated multi- 
ple frequencies on the COBE-XyMR experiment (three 
bands between 31 and 90 GHz), the current WMAP ex- 
periment (five bands between 23 and 94 GHz), and the 
future Planck experiment (nine bands between 30 and 
857 GHz). 

While the necessity of multi-frequency observa- 
tions has been recognized within the cosmological 
community for a long time, there has been un- 
certainty about how those observations should be 
utilized. Many different methods have been pro- 
posed, including the Maximum Entropy Method 
itBarreiro et al. 2004; Bennett et al. 2003b; Ho bson et alJ 
119981 iStolvarov et' alJ 12001 I2005D. the Inter n al Lin- 
ear Combina.tion method JBennett et al.l l2003bl: 
iTeemark et all 120031 lEriksen et al.ll2004aH . Wiener fil- 



tering l|Bouchet fc G isperi"1999': 'Tegmark fc EfstathioiJ 
19960: and the In dependent Component Analysis 
method llMa,ino et alJl2002t 120031 IDonzelli et a].ll200,qll . 
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Some of these methods (e.g., Baccigalupi et al. 2004 and 
StivoH et al. 2005) have been apphed to polarization 
data. These methods all have their origins in the field of 
computational image processing. 

In this paper, we advocate a more direct approach 
to the comp onent separat i on pro blem, following in the 
footsteps of iBrandt et alJ l)1994f) . We choose a phys- 
ically motivated parametric model for each significant 
signal component, and fit the free parameters by means 
of well-established algorithms, such as MCMC and non- 
linear searches. (For an alternative spectral matching al- 
gorithm, see Delabrouille, Cardoso & Patanchon 2003.) 
The advantages of this approach are many: it does not 
require priors (although it can benefit from them); no as- 
sumptions about the spatial structure of the foreground 
properties are imposed; the method scales proportion- 
ally to the number of pixels; the results may be verified 
by means of goodness-of-fit tests; and, most importantly, 
the method yields accurate uncertainties for all estimated 
quantities. 

We begin by developing a simple algorithm that is able 
to analyze real -world data in the presence of realistic 
noise. Whereas IBrandt et al.l l) 19941) relied exclusively on 
non-linear fitting, and was therefore quite unstable with 
respect to noise, we take advantage of more recent devel- 
opments (most importantly Markov Chain Monte Carlo) 
to make the algorithm both more robust, and also to 
produce accurate errors. We also suggest a procedure 
for propagating these results into final data products, 
namely the CMB power spectrum and cosmological pa- 
rameters, by means of a Gibbs sampling method. 

After establishing the algorithm, we apply it to realis- 
tic simulations corresponding to the future Planck and 
six-year WMAP data, taking into account the predicted 
noise distributions of each experiment. We attempt to 
model the foregrounds as accurately as possible, given 
our current understanding of the involved foregrounds. 

2. LIKELIHOOD FORMULATION OF THE COMPONENT 
SEPARATION PROBLEM 

Assume that the observed data take the form of a 
multi- frequency set of sky maps, d^, each of which may 
be written in the form 

Ay = As,, -I- n^, (1) 

where v identifies frequency bands and sky maps 1 
through N , is the true sky signal at the corresponding 
frequency band, is instrumental noise, and A denotes 
convolution with the instrument beam. 

We assume the noise component to be Gaussian dis- 
tributed with vanishing mean and variance a1{p), where 
p is the pixel number. Thus, the noise is uncorrelated 
both between pixels and between frequency channels, but 
spatial variations in the variance are allowed. 

Note that there is no frequency index on the beam 
operator in Equation ^ indicating that all channels are 
assumed to have the same beam response. For multi- 
resolution experiments, this implies that the sky maps 
must be smoothed to a common resolution prior to anal- 
ysis. Equation ^ is no longer strictly valid, since the 
noise is then correlated; however, in practice it works 
reasonably well to approximate the noise term as uncor- 
related between pixels, with rms levels determined by 
Monte Carlo simulations of processed noise. 



The signal may be decomposed into a sum of com- 
ponents, Su ~ X^i which the most important 
ones are the cosmological CMB signal and three galac- 
tic foregrounds; synchrotron, free-free and dust emission. 
Compact (unresolved) galactic and extra-galactic sources 
could also be included in the list, but these are more con- 
veniently detected by other methods, such as wavelets 
(e.g., Vielva et al. 2003). In this paper we consider com- 
pact source removal as a part of the pre-processing stage, 
and assume that resolved sources are either masked, me- 
dian filtered or fitted prior to analysis. Thus, we include 
only diffuse foregrounds in the following, but note that 
more work is needed on this issue. 

Assume also that the frequency spectrum of each sig- 
nal component may be parametrized by a small number 
of free parameters, for instance an amplitude and a spec- 
tral index S = A {y /vq)^ . Given a set of multi-frequency 
CMB sky maps as described above, and a parametric sig- 
nal model Sy{9) with free parameters 9, we now establish 

both a point estimate 9 for the free parameters, and the 
corresponding uncertainties. To do so, we use standard 
likelihood methods of parameter estimation. 

Since the noise is assumed to be Gaussian distributed 
and uncorrelated between pixels, the likelihood C reduces 
to x^, independent between pixels^'^. 



1 ^ 

2 ^ 



S.{9) 



1 



(2) 



The problem is thus reduced to mapping out a likelihood 
by some numerical technique, for instance grid computa- 
tion, MCMC, or non-linear searches. 

We illustrate the procedure in Figure ^ where the re- 
sults from an analysis of one arbitrarily chosen pixel are 
shown (see §|S1 for details). The observed data points 
are marked by black circles, and the fitted components 
are shown as smooth, colored curves. The dashed black 
curve shows the sum of all components. Four signal com- 
ponents are included in this model, CMB, synchrotron, 
free-free, and thermal dust emission. The parametric 
models for the foregrounds are perfect power laws for 
synchrotron and free-free emission, and a one-component 
model for dust (see Eauation ll2l) . 

Although we focus on temperature anisotropics in 
the present paper, the method can handle polarization 
anisotropics equally well. In that case, the takes the 
form 

N 

= E i^'^ - ^-(^)f (J, - S,{9)) , (3) 
1^=1 

where 

d.= \dQ\, S.{9) = SQ{9) (4) 
\d^J \S^{9)J 

are the Stoke's /, Q, and U parameters for the data and 
model respectively, and is the 3 x 3 (/, Q, U) noise 

In practice this is not strictly correct, since the CMB and 
foreground components are indeed correlated between pixels, and 
we will at some point work with smoothed sky maps, but it is a 
good approximat ion f or component separation purposes. Further, 
as described in S I5.6I at least spatial CMB correlations may be 
taken into account by Gibbs sampling. 
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Fig. 1. — Component separation using multi-frequency measurements (linear units in left panel, logarithmic units in right panel). Most 
signal components has a well-defined frequency spectrum that may be parametrized by one or a few parameters, and component separation 
may therefore be viewed as a standard parameter estimation problem. The example shown here is based on one single pixel in a simulated 
data set corresponding to the six-year WMAP and the Planck experiments, as discussed in §|H| The error bars on the data points are 
multiplied by a factor of 50 in order to make them visible on this scale. (Due to modeling errors, this particular fit has a of 44, and 
with five degrees of freedom, it is rejected at the 99.9999% confidence level.) 



correlation matrix for the pixel. (We still assume uncor- 
related noise between pixels and frequency bands, but 
not between the three Stoke's parameters for each indi- 
vidual pixel.) 

3. PARAMETRIC MODEL FITS 

In this section we introduce the parametric model that 
we fit to the data in the likelihood analysis, starting with 
a review of the currently favored parametric signal mod- 
els sj, for each signal component. We emphasize that the 
procedure as such is general, and any parametric model 
may be included in the analysis. Then we discuss how 
to take into account the effect of non-zero instrumen- 
tal bandwidths, and discuss implementation details that 
lead to more transparent computer code. Note that the 
parametric signal models we adopt for fitting the various 
components may differ from the parametric models that 
are used for modeling the components later in section lfTTl 
as will certainly be the case in dealing with real data. 

3.1. Signal components 

CMB anisotropies — The cosmological bac kground com- 
pone nt is, due to its black body nature ijMather et al.l 
most easily characterized by its thermodynamic 
temperature Tcmb, or equivalently, the anisotropy tem- 
perature ATniuR = Tc mr — Toj where Tq = 2.725 K 
ijFixsen fc Matheil l2002|) is the average CMB tempera- 
ture. However, as discussed later, the foreground com- 
ponents are more easily described in terms of antenna 
temperatures, and therefore we choose to convert the 
CMB signal accordingly. The CMB signal model then 
reads 

scmb(i^) = ArcMB _ -^y , (5) 

where x — hv/kTo, h is the Planck constant, and k is 
Boltzmann's constant. 



Synchrotron emission — Synchrotron emission from the 
Galaxy originates in relativistic cosmic ray (CR) elec- 
trons spiraling in the Galactic magnetic field. The mor- 
phology of the observed emission depends on the distri- 
bution of the relativistic electrons in the Galaxy, and the 
Galactic magnetic field structure. In the Galactic plane, 
the latter exhibits a large-scale ordering with the field 
parallel to the spiral arms (the regular component). Su- 
perimposed on this is real small-scale structure (the irreg- 
ular component) which shows variations between the arm 
and inter-arm regions and with gas phase. The regular 
and irregular components seem to be of comparable mag- 
nitude. At high latitudes, there is a contribution from the 
Galactic halo, and specific nearby structures (e.g., the 
North Polar Spur). Variations in the frequency spectral 
index of the synchrotron continuum emission arise from 
variations in the CR electron energy spectrum, which has 
a range of distributions depending on age and the envi- 
ronment of origin (e.g., supernova explosions or diffuse 
shocks in the interstellar medium). 

The synchrotron emission may be accurately modeled 
by means of a simple power law over a considerable range 
of frequencies, 

s,{i.)^aJ—Y\ (6) 

Here is the synchrotron amplitude^'* (measured in an- 
tenna temperature /iK) at some reference frequency z^o,s, 
and /3s is the synchrotron spectral index. Since the spec- 
tral index varies with both frequency and position on the 
sky, at least two free parameters are required to describe 
the synchrotron emission properly in a given direction. 

fLawson ct al. ( 1987) studied the spectral index vari- 
ation based on low- frequency radio surveys, and found 
that the brighter regions away from the Galactic plane 

Note that italic font A is used for component amplitudes, 
while bold font A is used for beam convolution. 
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have typical valu es of /3 at 100 and 800 MHz of 2.55 and 
2.8, respectively. iReich fc ReichI l)1988j) used radio con- 
tinuum surveys of the Northern sky at 408 and 1420 MHz 
to demonstrate a range of spectral index values between 
2.3 and 3.0, with a typical dispersion A/3 = ±0.15. The 
steepest spectra were observed towards the North Polar 
Spur, and there was a flattening in spectral index to- 
wards higher latitudes in the Galactic anti-centre direc- 
tio n. Such behavior h as been confirmed over the full sky 
bv IReich et al.l l)200 3l . who find that spectral fiattening 
is particularly pronounced in the Southern hemisphere 

At higher frequencies, the brightness temperature 
spectral index is expected to steepen by ^ 0.5 due to elec- 
tron e nergy losses (Platania et al.llT99^ . iBandav et al.l 
(|2P03') derived a mean spectral index between 408 MHz 
and 19.2 GHz from the Cottingham ( 1987) survey and be- 
tween 31.5, 53, and 90 GHz from the COBE-BMR data. 
The steep spectral index of ~ 3.1 for Ga lactic latitudes 
| 6| > 15 ° is consistent with expectations. iBennett et al.l 
l)2003tf > claim that the spectral break occurs near the K- 
band. Spectral indices above 10 GHz are likely between 
2.7 to 3.2. 

Free-free emission — Free- free emission is the brems- 
strahlung radiation resulting from the Coulomb interac- 
tion between free electrons and ions in the Galaxy. Free- 
free emission is weaker than synchrotron emission at fre- 
quencies below ^ 1 GHz over most regions of the sky, 
with exceptions in the Galactic ridge and the brighter 
featu res of the local Gould Belt system ( Dickinson ct al. 
I2003j) . Free-free almost never dominates at high latitudes 
in any frequency band, and is therefore difficult both to 
observe and to simulate. 

Usin g the relations summarized in 'Dickinso n et al.l 
l)2003tl . it can be shown that the brightness temperature 
of the free-free emission is described by the relation 



Tff,b oc iy-'^T-°-^ (ln[0.04995z/-i] + l.SlnTo 



(7) 



w here Tp, is the electr on temperature. 

IShaver et al.l ljl983f) used radio recombination lines to 
show that the electron temperature of Hii regions at the 
galacto centric radius of the Su n (i?o = 8.5 kpc) is 7200 ± 
1200 K. IPaladini et alJ l|2005j) found similar results from 
a larger sample containing many weaker sources. At high 
galactic latitudes, the ionized hydrogen typically will be 
within 1 kpc of the Sun (,Dickinson et„aL. .2QQ3,) , and 
we expect the electron temperature to be in the range 
Te = 7000-8000 K, although it is possible that the diffuse 
emission at a given Galactocentric distance may differ 
from that of the higher density Hii regions on the plane. 

The frequency dependence is therefore well con- 
strained, with an effective spectral index of Ps = —2.14 
at the frequencies of interest, the range over 10-100 GHz 
being of order —2.1 to —2.2 and steepening still further 
t o —2 .3 at hundreds of GHz. 

iFinkbciner ( 2004a) has analyzed the WMAP data and 
found a significant component with a free-free like spec- 
trum within 30° of the Galactic center. The component 
is uncorrelated with Ha emission and may be indicative 
of hot 10^ K) gas. We do not attempt to include such 
a component here. 

In this paper, we model free-free emission as a simple 
power law with a fixed spectral index, 



Sff(j/) = As [ 



(8) 



Thus only one parameter is required for free-free emis- 
sion. Future experiments may need to estimate the elec- 
tron temperature directly from the data. In that case 
Equation [7| should be used directly, at the cost of intro- 
ducing one extra free parameter into the fit. 

Thermal dust emission — The thermal dust emission that 
contributes to the frequencies of interest for CMB anal- 
ysis arises from grains large enough to be in thermal 
equilibrium with the interstellar radiation field, and is 
known from analysis of the IRAS and COBE-BIRBE 
data to peak at a wavelength of approximately 140 /im. 
At higher frequencies, there is a contribution from the 
optically active modes of PAH molecules, but these are 
not of interest here. 

Currently preferre d dust emission models 
l|Finkbeiner et al.ll9 99^ extrapolate from high-frequency 
COBE-FIRAS and -DIRBE observations to CMB fre- 
quencies using combinations of modified blackbody fits 
and accounting for dust temperature variations. Such 
fits approximate the integrated contributions to the 
emission from multiple components of dust, i.e., with 
different grain properties (chemical composition and 
size) and equilibrium temperatures. The best-fit model 
(model 8 of Finkbeiner et al. 1999) assumes two main 
components: 



(9) 



Here F represents the combined COBE-DIRBE and 
IRAS template (Schlcgel ct al. 1998), i^o,d = 3000GHz, 
and l3d{i^) is dependent on the frequency as discussed 
above, 

/?d(.) = , (10) 



d{y) 



91 

92 



V 



3 + Ql 



/2 



fO,d 



3+Q2 



(11) 



e '--^2 - 1 



with best-fit parameters /i = 0.0363, <7i/(j2 = 13, ai — 
1.67, = 2.70, T i = 9.4 K, T2 = 16.2 K, and /2 = I-/1 
iFinkbeiner eTall (|T999') . 

In principle, these equations may serve as our paramet- 
ric model for fitting the dust emission spectrum. How- 
ever, few (current or future) CMB experiments have suf- 
ficient power to constrain six parameters for dust alone, 
and simplifications are therefore unavoidable. Rather 
than fitting the full form as g iven above, we therefor e 
choose the simpler "model 3" of lFinkbeiner et all l)1999j) . 
setting /i = 1 and Ti = 18.1 K, but letting ai vary freely. 
Equation El may then be simplified to 



Sd{v) = ^d- 



e '=^1 



V 

Vd.Q 



(12) 



where Ad is the thermal dust amplitude at a reference 
frequency I'd.O) s-^d a /3d is a free parameter. Thus, the 
fitted model is a power law modulated by a slowly de- 
creasing function of order unity over the frequencies of 
interest. 



5 



Such a parameterization does not allow for the spec- 
tral break that does exist in the COBE-FIRAS at ap- 
proximately 500 GHz, and which may reflect either the 
emissivity of different grain components, frequency de- 
pendence of the emissivity of the dominant grain compo- 
nent, or possibly a pop ulation of cold dus t grains mixed 
with the warmer dust ijReach et alJll995|l . However, at 
500 GHz the signal of dust anisotropics is so strong com- 
pared with CMB anisotropics that it is useless for pur- 
poses of component separation. 

Anomalous dust emission — Cross-correlation of the 
COBE-BMR data with the DIRBE map of thermal dust 
emission at 140 /im in Kogut ct al. (199fD revealed an 
anomalous component with rising spectrum from 53 to 
31.5GHz. Bandavct al. (2003), again using the DMR 
data together with a survey at 19.2 GHz, and indepen- 
dently the WMAP team ( Bennett et al. 2003b), sug- 
gested that this component was well-described by a 
power-law spectrum with index —2.5 for frequencies in 
the range ~20-60GHz. The latter proposed that the 
emission originates in star-forming regions close to the 
Galactic plane. However, the favored model to explain 
this anomalous dust-correlated component is in terms of 
t he rotational emission from very small grains. 

iDraine fc LazarianI 1)19981) have developed a three com- 
ponent model of this 'spinning dust' which contains 
contributions from the three phases of the interstellar 
medium — the Cold Neutral Medium, the Warm Neutral 
Medium, and the Warm Ionised Medium. The charac- 
teristic spectral behavior of the model includes a rising 
spectrum up to a turn-over in the range 10-20 GHz, then 
a rapidly falling spectrum which can be characterized by 
an effective spectral inde x in excess of 3 bey ond 30 GHz. 

Recent observations bv lFinkbeinerl l|2004bj) of dust cor- 
related emission outside of Hli regions between 8 and 
14 GHz shows a rising spectral slope and amplitude far 
exceeding that as sociated with therm al dust emission. 
More importantlv. lWatson et ahl l)2005)) show that obser- 
vations of the Perseus molecular cloud made between 11 
and 17 GHz and augmented with the WMAP data can be 
adequately fitted by a spinning dust model. Neverthe- 
less, although the case is compelling for describing the 
anomalous emission by such models, the detailed mor- 
phology of the emission remains uncertain and no un- 
ambiguous template to trace it exists. It is also possible 
that the emission may be confined to specific clouds at 
relatively low Galactic latitude, leading to a more patchy 
distribution than for the diffuse thermal dust contribu- 
tion with some additional sensitivity to environment. In- 
deed, iLagachel 1)2003)) presents evidence that the excess 
emission is associated with small transiently heated dust 
particles, which may be destroyed under certain physical 
conditions. Given these uncertainties, we do not include 
this foreground component in our studies. 

3.2. Non-zero bandwidths 

The previous sections describe the basic behavior of 
each signal component as a function of well-defined fre- 
quencies. However, real experiments integrate over a 
range of frequencies (a "frequency band"), typically with 
unequal weights, and the observed signal strength does 
not equal that given by the central frequency alone. 

We take into account this effect through the concept 



of an "effective frequency" z^cff, defined by 

SM = [ f{iy)S,diy, (13) 



where S^, is the frequency spectrum of the signal, and 
f{v) is the frequency response profile of the detector. 
Thus, the spectrum at the effective frequency equals the 
average over the frequency band. The advantage of doing 
this is simply that computationally expensive integrals 
are replaced by single point computations. 

In this paper we assume for simplicity that all fre- 
quency response functions correspond to flat bandpass 
filters with sharp frequency cutoffs at i>a and i>h. For 
simple power law models, such as those of synchrotron 
and free-free emission as described above, the effective 
frequency of a signal component with spectral index P is 
then given by 



1 



,/9+l 



,/3+l 



1//3 



(3+1 Vb-Va 



(14) 



For more complicated spectra. Equation 1131 must be 
solved numerically. Fortunately, it is straightforward to 
pre-compute a grid of the effective frequencies prior to 
the full analysis, since they only depend on the frequency 
scalings and not the component amplitudes, and compu- 
tational speed is not compromised. 

3.3. Implementation details 

To simplify the computer code, it is convenient to in- 
troduce some general notation. For instance, if we can 
write all signal models in a common form, we do not have 
to consider a list of special cases, but rather handle all 
cases with the same code. 

Indeed, all frequency spectra discussed above may be 
written in a common form, namely that of a power law 
modulated by an arbitrary frequency-dependent func- 
tion, 



(15) 



Here Ai {p) and (3i {p) are the "free" amplitude and index 
parameters for component i in each pixel p, and Q^^^ff. is 
an arbitrary function only dependent on frequency. 

Specifically, synchrotron and free-free emission are in- 
cluded simply by setting Ci^i, — 1, while the CMB compo- 
nent is defined by /3cmb = and Ccmh,u — x^e^/ {e^ — 1)^, 
as discussed above. For dust, Cd i, is given by equation 

112 

Even anomalous dust could be included within this no- 
tation. One option is simply to tailor the correction fac- 
tors c,:.,^ to match the pre dicted spinning dust spectrum 
l)Draine fc LazariarJl998D . fix the corresponding spectral 
index at zero, and then fit for the amplitude only. An- 
other is to merge the spectrum with that of the thermal 
dust emission, and thereby enforce identical spatial tem- 
plates. 

4. SUMMARY OF PREVIOUS RESULTS 

This paper may b e seen as a na t ural c ontinuation of 
the work started by iBrandt et al.l 1)19941) , who consid- 
ered how well future experiments could reconstruct the 



6 



CMB signal in presence of noise and foregrounds. Their 
approach, parameter estimation, was the same as ours; 
however, they relied solely on maximum likelihood es- 
timation (i.e., non-linear fitting), and their results were 
therefore less stable with respect to noise than the ones 
we present here, as will be seen below. Nevertheless, sev- 
eral of the conclusions drawn bv lBrandt et al.l l)1994|) are 
still valid for our work, and well worth repeating: 

1. The number of frequencies must equal or exceed 
the number of fitted parameters, else the problem 
is mathematically degenerate. This is obvious, but 
not trivial: no experiment to date has had the min- 
imum number of frequencies required to separate 
CMB, synchrotron, free-free, and dust fluctuations, 
even in their simplest form. 

2. One should attempt to reduce the number of free 
parameters in the problem, as this gives greater sta- 
bility with respect to noise. Seemingly gross simpli- 
fications, such as approximating both synchrotron 
and free-free emission by a single power law, can 
often yield improvements in the reconstruction. 

3. It is usually advantageous to fit spectral parameters 
to reduced-resolution and low-noise data, and then 
solve for the amplitudes in the full-resolution data, 
fixing the indices at the smoothed values. 

4. Due to the similarity between the synchrotron 
and free-free emission, better results are obtained 
whenever the latter is not a significant contami- 
nant. Thus, if the free-free contamination could be 
constrained by radiation physics knowledge, it is 
well worth trying. 

While working on the present analysis, we have repro- 
duced all of these conclusions, and most of them have 
been taken into account when establishing the prescrip- 
tion described below. However, we will not elaborate 
further on thes e issues here, but ra ther refer the inter- 
ested reader to iBrandt et al.l l|1994|) , and come back the 
above points as they are needed in the analysis. 

5. METHOD 

In this section we propose an algorithm for solving the 
parameter estimation problem with sufficient speed and 
accuracy to be useful for practical analysis of current 
and future data. Each step of the algorithm consists of 
well-established methods, and the approach should seem 
quite familiar. 

5.1. Overview 

The goal is to establish both a point estimate of all 
interesting parameters, and their uncertainties. Our pre- 
scription for doing so is as follows: 

1. If required by noise levels or computational re- 
sources, downgrade all sky maps both in pixel and 
beam resolution. 

2. For each low-resolution pixel, 

(a) choose a parametric model, 

(b) solve for all parameters jointly by MCMC, 



(c) estimate non-linear parameters and corre- 
sponding uncertainties by marginalizing over 
all other parameters, 

(d) find the goodness-of-fit in terms of the ■ 

3. For each high-resolution pixel within a low- 
resolution pixel, 

(a) either (fast but approximate analytical ap- 
proach) 

i. fix the non-linear parameters at the low- 
resolution pixel values, 

ii. find maximum-likelihood estimates for all 
linear parameters (i.e., component ampli- 
tudes) by solving a linear equation, 

iii. find corresponding uncertainties by ana- 
lytic error propagation formula, 

iv. estimate the goodness-of-fit in terms of 
the again. 

(b) or (exact but expensive Gibbs sampling ap- 
proach) 

i. use the low-resolution MCMC distribu- 
tions to sample non-linear parameters nu- 
merically, 

ii. given a set of non-linear parameters, sam- 
ple amplitudes from their corresponding 
Gaussian distribution, 

iii. given foregrounds, sample CMB sky map 
and power spectrum. 

The route outlined in step 3b) holds promise of a com- 
plete solution to the foreground problem in a CMB con- 
text, since, if successful, the foreground uncertainties 
are propagated all the way from noisy observations to 
the CMB power spectrum and cosmological parameters. 
However, in this paper we only present the basic ideas, 
and leave the details for a future more comprehensive 
study. All high-resolution results presented in the fol- 
lowing are thus based on the analytical approach. 

5.2. Non-linear parameters and large-scale smoothing 

One of the main themes of lBrandt et alJ l|1994f) was the 

instability of a non-linear fit with respect to noise. This 
is not hard to understand. In estimating multiple pa- 
rameters from a limited number of frequency bands, the 
maximum likelihood point may easily slide along some 
degeneracy ridge on the likelihood surface in the presence 
of realistic noise. For all currently planned experiments, 
additional degree-scale smoothing is a requirement in or- 
der to reach reconstruction errors at the /iK level. 

Another and more practical issue is the fact that our 
main algorithm relies on MCMC analysis of each indi- 
vidual pixel. This takes on the order of 100 seconds per 
pixel. Even though the algorithm parallelizes trivially 
because the pixels are analyzed individually, a complete 
analysis at full Planck resolution 50 million pixels) 
would be unfeasible. 

It is important to realize, however, that a full MCMC 
analysis is required only for estimating non-linear pa- 
rameters, such as spectral indices or dust temperatures. 
If all parameters in the problem are linear (i.e., compo- 
nent amplitudes), an analytical computation is equally 
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good. We therefore compute the complete probabil- 
ity distributions from reduced-resolution maps, fix the 
high-resolution non-linear parameters at the correspond- 
ing low-resolution values, and then solve for the high- 
resolution component amplitudes with alternative meth- 
ods, for instance analytically or by Gibbs sampling. 

When adequate data on foregrounds are in hand, the 
validity of this smoothing can be tested. If preliminary 
indications turn out to be true, and spectral indices vary 
more slowly on the sky than amplitudes, the smoothing 
process will not lead to significant loss of information. In 
any case, the smoothing scale can be optimized within 
the bounds of computational resources. 

5.3. Model selection 

A second main theme of iBrandt et alJ l)1994f) was the 
importance of model selection. They clearly demon- 
strated that a large number of parameters does not nec- 
essarily yield a better CMB reconstruction. Quite the 
contrary, too many parameters often yield unphysical re- 
sults. In general, one should never fit more parameters 
than required by the data. 

In principle, it would be useful to have an automated 
prescription to identify the optimal model for a given 
pixel. To some extent, such a procedure is provided by 
means of the so-called information criterio n (IC). an ide a 
that was introduced to CMB analysis bv iLiddla l)2004|) . 
The fundamental idea in this approach is not to maximize 
the likelihood alone, but rather minimize the IC, defined 
as follows, 

lC^-2\nC + ak. (16) 

Here k is the number of parameters in the fit, and a is 
a penalty factor. (Two useful choices for a are aA = 2, 
the Akaike information criterion, and ckB — InA^, the 
Bayesian information criterion, N being the number of 
data points.) Within this framework, a new parameter 
must prove its usefulness by returning a significant im- 
provement in the fit to be included in the model. 

We implemented this approach in our codes, and ob- 
tained reasonable results. However, the model sky map 
had a clear tendency to be patchy, and not necessar- 
ily well correlated with physical structures. In the cur- 
rent implementation, we therefore only use the infor- 
mation criterion approach to inform our model choices, 
and tailor the model map manually according to known 
structures. For instance, in the example given in §|H1 
we use a full four-component model (CMB, synchrotron, 
free- free, and dust) insid e an expanded KpO galactic 
cut l)Bennett et al.l[2003bj) . as well as in a few selected 
patches (e.g., the LMC), but we ignore free-free other- 
wise. However, we expect that the information crite- 
rion approach may be developed further, and should be 
a valuable tool for future experiments. 

Realistically, model selection is likely to be an iterative 
process, as will be demonstrated in the worked example 
of §0 Typically, one constructs an initial physically mo- 
tivated model map, and performs the analysis with that 
model set. Based on the results, one then evaluates the 
goodness-of-fit for each pixel (i.e., the as defined by 
equation 121, and compares the results with a simple 
distribution with the appropriate number of degrees of 
freedom (number of frequencies minus the number of free 
parameters). If the agreement is poor, the model set can 



be modified, or offending pixels can be rejected from fur- 
ther analysis. 

Using such methods, model selection is likely to be a 
relatively straightforward (although somewhat tedious) 
process for future experiments. However, in the present 
paper we are satisfied with a very simple choice of mod- 
els, based on established sky masks. 

5.4. Parameter estimation by MCMC 

By now, parameter estimation by Markov Chain Monte 
Carlo is a well established technique within the CMB 
community, with its most visible application being esti- 
mation of cosmological parameters from the CMB power 
spectrum. To our knowledge, it has not yet been ap- 
plied to component separation, and we therefore briefly 
describe the algorithm here. For more details , we r efer 
the interested rea der to, e.g., iLewis fc Bridlel l)2002|) or 
IVerde et aki ll2nnl. 

5.4.1. Algorithm 

Suppose we want to estimate a set of parameters and 
corresponding uncertainties from a set of observed data, 
and that we know how to compute the likelihood given an 
arbitrary combination of parameters. The Markov Chain 
Monte Carlo algorithm is then given by the following 
simple steps: 

1. Choose any initial point in parameter space, Oq, 
and compute the corresponding likelihood, Cq — 

2. Define a stochastic function / that, given param- 
eters 9i, returns a new set of parameters fl^+i = 

f{0^)■ 

3. Compute 9i+i given by / and the corresponding 
likehhood, £i+i — £{9i+i). 

4. Set 9i+i — 9i (i.e., reject the proposal) with prob- 
ability p = 1 — min(£i+i/£i, 1). 

5. Go to step no. 3, and iterate as long as necessary. 

This procedure returns a chain of parameter samples 
9i, i = 1, . . . , iVsampies, and their multi-dimensional his- 
togram equals the likelihood in the limit of an infinite 
number of samples. 

Some intuition for the process may be gained by notic- 
ing the form of the "jump probability" given in step 4 of 
the algorithm: If the likelihood of the proposed point is 
larger than that of the old point, we never reject the 
proposed point; we always move towards more likely 
solutions when proposed. However, if the likelihood is 
smaller, we still accept the proposed point with proba- 
bility p = Ci+i/Ci. This guarantees that we spend most 
of the time around the peak position, but still explore 
less likely points. Indeed, it may be proven that the time 
spent at a given point in parameter space is proportional 
to the likelihood itself. 

5.4.2. Automated MCMC in practice 

In practice, there are several problems connected to 
MCMC parameter estimation; usually, most of these may 
be identified (and often solved) by simple visual inspec- 
tion of the Markov chains. However, since we want to 
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analyze many thousands of independent pixels, finding 
automated and yet robust solutions to the same prob- 
lems is of critical importance; a solution that works in 
99% of all cases is not good enough. 

Burn-in — Although it is true that the initial guess may 
be chosen arbitrarily, and that the chain eventually will 
burn in to the right solution, it is difficult to construct 
a truly reliable automated prescription for when burn- 
in has occurred. Also, since computational speed is of 
critical importance, it is not acceptable to spend a long 
time in the burn-in phase. For both of these reasons, we 
choose to initialize the chains at the maximum-likelihood 
point, which we find using a standard non- linear fitting 
algorithm. (We have found that a sequential quadratic 
programming (SQP) method works very well for this 
task.) 

Proposal function — In step 2 of the algorithm, we must 
establish a proposal function /. A simple example of 
such a function would be 6l_^_i — 9^ + SB^rj^ , where j is a 
parameter index, 69^ is a pre-defined rms step size, and 
is a Gaussian stochastic variate of zero mean and unit 
variance. However, since most parameters of interest are 
usually strongly correlated, this choice is quite inefficient. 

Our current best solution is to run a preliminary chain 
(using the univariate Gaussian proposal function de- 
scribed above, with manually set step sizes), and com- 
pute the covariance matrix Cij = (59^59^^ of the result- 
ing samples. We then Cholesky-decompose this matrix, 
C = LL', and define our new proposal function to be 
Bi+i = 0i-|-aL?7, where rj now is a vector of Gaussian vari- 
ates and a is an overall scale factor, typically initialized 
at ~ 0.3. This ensures that the proposed samples have 
approximately the correct covariance structure, and the 
overall sampling efficiency is thereby greatly improved. 

To avoid too large or too small step sizes, we also im- 
pose the requirements that the acceptance ratio (the ra- 
tio of accepted to rejected proposals) must be higher than 
5%, and lower than 80%. If one of these two criteria is 
violated, we divide or multiply a by 2, and re-start the 
MCMC analysis. 

Convergence — Finally, we must decide when a sufficient 
number of samples has been accumulated. No general 
solutions are available. We have adapted a good work- 
ing solution proposed by Gclman & Rubin (1992), as fol- 
lows. Run m independent MCMC chains in parallel for 
the same pixel, each producing n samples. Compute the 
following quantities, 

»' = ;;;^;^EE(''; -»0 • m 

V = (l--]w+-B, (19) 
\ n J n 

R-^, (20) 

where 9 is the average over all m-n samples, and 9^ is the 
average of the samples within chain number j. W esti- 
mates the variance within each chain individually, while 



B estimates the variance between the chains. When the 
chains have converged properly, V and W should be iden- 
ti cal, a nd R should be close to unity. 

iGelm an & Rubin ( 1992) make the general recommen- 
dations that the initial points for the m chains should 
be over-dispersed relative to the true distribution, and 
that the chains should be run until R < 1.2. However, 
as discussed above, we initialize the chains at the max- 
imum likelihood value in order to avoid burn-in prob- 
lems, and thus the first point is certainly not fulfilled 
in our approach. Consequently, the numerical value of 
the convergence criterion they give does not apply to our 
prescription. 

To remedy this situation, we impose two alternative 
criteria. First, we require that the chains run for a min- 
imum number of samples (typically on the order of 10^, 
but only storing every, say, 500th sample, to reduce cor- 
relations). Second, we require that the largest value of 
R, individually computed for all parameters included in 
the model, must be smaller than 1.01. In most cases, we 
find that the latter criterion is fulfilled long before the 
former, indicating that the overall CPU time may be de- 
creased somewhat. However, since we cannot inspect the 
chains manually for more than a few pixels, we consider 
safety to be more important than speed, and adopt a 
very conservative approach. With the criteria discussed 
here, we have found excellent convergence in all cases we 
have inspected, and the Monte Carlo error (the error on 
the error due to a finite number of samples) is typically 
less than 1%. 

5.4.3. Point estimators and uncertainties 

The MCMC algorithm provides us with a large num- 
ber of multi-variate samples drawn from the likelihood, 
and these may be used to form a great variety of use- 
ful statistics. Here we focus on the uni- variate distri- 
butions for each parameter, marginalizing over all oth- 
ers. (Marginalization with MCMC samples is straight- 
forward: simply disregard the "uninteresting" parame- 
ters, and make a histogram of the "interesting" parame- 
ter sample values. Priors are discussed briefiy in the next 
section.) 

Our point estimate for each parameter value is then 
defined to be the mean of the MCMC samples, with 
an uncertainty given by variance of the samples. As 
shown later in practical examples, this Gaussian approx- 
imation is quite good for both the CMB temperature 
and the spectral indices, while the distributions for the 
synchrotron and free-free amplitudes tend to be non- 
Gaussian due to a combination of strong internal cor- 
relations and a positivity prior. 

One of the most importa nt differences between our ap- 
proach and that taken by iBrandt et al.l l)1994f) is that 
they chose the maximum likelihood value as their point 
estimate, whereas we choose the mean. This makes our 
estimate considerably more stable with respect to noise, 
since it takes more to shift the entire likelihood volume 
than to change its shape. Therefore, not only does our 
method yield accurate error bars on all relevant quanti- 
ties, but also the point estimates are more reliable. 

5.4.4. Priors 
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Bayes' theorem states that the posterior distribution 
P{6\d) is given by 



p{0\d) oc p{d\e) p{e) = c{e) p{e), 



(21) 



where £{6) — P{d\6) is the hkehhood and P{0) is a prior. 
ImpUcit in the above discussion, we have adopted the 
simplest possible choice for this paper, namely uniform 
priors between two (not necessarily finite) limits. Specifi- 
cally, we impose no constraints on the CMB temperature, 
and only a positivity constraint {A > 0) on the fore- 
ground amplitudes. For the spectral indices, our priors 
are sufficiently generously chosen to never exclude phys- 
ically realizable values (—4 < f3s < —2.2 for synchrotron, 
and 1 < /3d < 3 for thermal dust). However, anticipating 
a future WMAP-only analysis, priors may then become 
considerably more important due to the high noise levels 
of the experiment. 

5.5. Analytic estimation of linear parameters 

The low-resolution MCMC analysis is by far the most 
expensive step in our algorithm. Once this step has 
been completed, estimation of the component amplitudes 
may be done both accurately and efficiently by solving 
sets of linear equations, with negligible computational 
resources. 

Recall that our likelihood function takes the form of a 
standard x^i 

2 



X 



dh 



fc = l 



Ncomp 

E 

2=1 



Ai Cik ( — 



13, 



(22) 



Since the only free parameters in the problem are now the 
component amplitudes A^, and the noise is assumed to 
be Gaussian, this is simply a multi-variate Gaussian dis- 
tribution. Therefore the mean of the distribution equals 
the maximum likelihood value, and may be determined 
simply by equating the derivatives of the with respect 
to the parameters to zero. In a matrix form, this reads 







da 

where a' — {Ai, . . . , An^^ 



a = M 



E 

k=l 



Cik Cjk 



Vjk_ 



and 



E 



Via 



(23) 



(24) 



(25) 



While the above formulae yield excellent results for the 
high-resolution parameter point estimates (as long as all 
amplitudes are non- negative) , reliable error estimation 
within this framework is complicated. One problem is 
that we need to propagate the errors in the spectral in- 
dices into the final data products, taking into account the 
strong correlations between the errors. A second problem 
is introduced by the positivity prior on the foreground 
amplitudes, which leads to strongly non-Gaussian dis- 
tributions when active. Nevertheless, for well-behaved 
pixels (i.e., those with clear detection of all components 
individually) a rough approximation may be established 
by means of the usual error propagation formula. 



Suppose we are interested in a quantity z = f{x, y, . 
that depends on a set of measured quantities x,y,. 
each with independent and Gaussian errors Ax, Ay, . 
In this case, the uncertainty Az is given by 



Az' 



dl 

dx 



Ax' 



dl 
dy 



Ay2 + 



(26) 



This may be applied to our case by making the iden- 
tification Ai — /(x) = Y^^^r^ ^ij^dj. The uncertain 
quantities are both the observed data and the non-linear 
parameters, x* = {d^,(3Y. 

To compute the uncertainties, we need the partial 
derivatives, which by equation 1231 read 

|^ = M-^M-M + M-|^. (27) 
ax ax ax 

The derivatives 9M/9x and dd/dx are obtained from 
equations |^ and 123 

Great care must be taken when applying this method 
to the high-resolution data — it is strictly valid only under 
the assumptions that the uncertain quantities are both 
Gaussian distributed and internally independent, neither 
of which is true for our problem. Nevertheless, while the 
formal requirements are not strictly fulfilled, the approx- 
imation may still be useful for establishing the order of 
magnitude of the uncertainties. 

First consider the assumption of independence. The 
data points are for all practical purposes uncorrelated 
with the non-linear parameters, since the latter are esti- 
mated from heavily smoothed maps, and therefore aver- 
aged over many pixels. Additionally, if the model only 
contains one spectral index, or two weakly coupled in- 
dices (e.g., synchrotron and dust indices, but not syn- 
chrotron and free- free), the independence assumption is 
thus nearly fulfilled. 

The assumption of Gaussian distributions is generally 
more problematic. On the one hand, if a particular fore- 
ground component is near or below the detection limit, 
the actual distribution is sharply cut off for negative val- 
ues. On the other hand, if two components have simi- 
lar spectral dependencies (like synchrotron and free-free 
emission), the corresponding parameter probability dis- 
tributions are usually strongly skewed. In either case, 
the rms is not an accurate estimator of the uncertainty. 

5.6. Propagation of errors to the CMB power spectrum 
and cosmological parameters 

The analytical approach described in the previous sec- 
tion yields good point estimates for the desired param- 
eters, but only approximate uncertainties. Further, it is 
not straightforward to propagate the errors further into 
higher-level data products such as the CMB power spec- 
trum or cosmological parameters. 

A much more powerful solution may be devised by 
combining the methods described in the present pa- 
per w i th the Gibbs s ampling ap proach o f Jewell ct ajj 
(j200^ . rWa^delt et alJp304t) and lEriksen et al. (2004^)- 
Whereas most other techniques only provide the user 
with a very simple description of the power spectrum 
probability distribution (e.g., a maximum-likelihood es- 
timate and a Fisher matrix), the Gibbs sampling ap- 
proach yields the complete multi-variate probability den- 
sity P{Ci\d), Ce being the CMB power spectrum and 
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TABLE 1 

Frequency band specifications 



Experiment 


Center frequency 
(GHz) 


Bandwidth 
(GHz) 


Beam FWHM 
(arcmin) 


RMS per 7' pixel 

(mK) 


RMS per 60' beam 

(mK) 


WMAP 


23 


5 


52.8 


50 


± 7 


7.6 ± 1.1 


LFI 


30 


6 


33.0 


33 


± 7 


3.0 ± 0.6 


WMAP 


33 


8 


39.6 


51 


± 7 


5.0 ± 0.7 


WMAP 


41 


11 


30.6 


49 


± 8 


4.3 ± 0.7 


LFI 


44 


8 


24.0 


33 


± 6 


2.7 ± 0.5 


WMAP 


61 


16 


21.0 


60 


± 9 


4.8 ± 0.7 


LFI 


70 


14 


14.0 


31 


± 6 


2.4 ± 0.5 


WMAP 


94 


24 


13.2 


73 


± 11 


5.6 ± 0.8 


HFI 


100 


33 


9.5 


14 


± 3 


1.1 ± 0.2 


HFI 


143 


48 


7.1 


8 


± 1 


0.6 ± 0.1 


HFI 


217 


72 


5.0 


11 


± 2 


0.9 ± 0.2 



Note. — Specifications for each frequency band used in the simulation. All beams are assumed to 
be Gaussian. The RMS values for 7' pixels are computed taking into account the scanning strategy of 
each detector, but neglecting noise correlations. The RMS values per 60' beam are estimated from 1000 
Monte Carlo simulations by drawing Gaussian random numbers corresponding to the RMS level of 7' 
pixels, deconvolving the instrument beam, and finally convolving with a 60' FWHM Gaussian beam. The 
values also take into account reduced pixel resolution, from 7' to 14' pixels. The Planck RMS values are 
requirement levels, not goals. 



d the data. Further, it is very straightforward to in- 
troduce new sources of uncertainty into the framework, 
and such uncertainties are then seamlessly propagated 
through to the final data products. We outline here how 
foreground uncertainties may be propagated to the CMB 
power spectrum and cosmological parameters, but leave 
the details for a future publication. 

The Gibbs sampling approach is similar in philosophy 
to the MCMC method that is used extensively in this 
paper: the target density is established by drawing sam- 
ples from it. In our case, we are interested in the joint 
probability distribution P(C£, scmb, Ss, sg, Sd|d), where 
Si are the four interesting signal components discussed 
earlier, CMB, synchrotron, free-free, and thermal dust 
emission. While it is difficult to sample from this distri- 
bution directly, the Gibbs sampling algorithm provides 
a neat solution. Suppose we want to draw samples from 
a joint distribution P{x,y), but only know how to sam- 
ple from the conditional densities P{x\y) and P{y\x). In 
that case, the theory of Gibbs sampling says that samples 
(x, y) can be drawn by iterating the following sampling 
equations, 

x^+i ^ P{x\y'), (28) 
^ P{yW+'). (29) 

The symbol ' indicates that a random number is 
drawn from the distribution on the right hand side. After 
some burn-in period, the samples will converge to being 
drawn from the required joint distribution. 

Suppose now that we want to analyze a data set that 
for simplicity only includes CMB and synchrotron emis- 
sion, the latter being parameterized by an amplitude As 
and a spectral index /3s for each pixel. Suppose further 
that we already have run an MCMC analysis for each 
pixel as described earlier, and have access to the cor- 
responding probability distributions. In that case, the 
Gibbs sampling algorithm may be applied by means of 



the following sampling chain: 

/3:+i^P(/3.|Q,A:,s^CMB,d), (30) 
A:+i^P(^|Q,/3:+i,sbMB,d), (31) 
ScMB - P{scMB\ClAl+\Pl+\d), (32) 
q+i ^ P(Q|sj+jB)- (33) 
(The CMB power spectrum only depends on the CMB 
signal, not the foregrounds, and therefore the other com- 
ponents are omitted from the right-hand side in the last 
equation.) The first two rows are to be performed for 
each pixel individually, while the last two rows are per- 
formed in harmonic space, reflecting the intuitively pleas- 
ing idea that foregrounds should be handled in pixel 
space, while CMB fluctuations are better handled in har- 
monic space. 

To perform the analysis as described above, we have 
to be able to sample from all involved conditional dis- 
tributions. Sampling the CMB signal and power spec - 
trum parts is detailed by, e.g., lEriksen et a D ll2004bD . 
Sampling the foreground amplitudes given the spectral 
indices is straightforward, since the corresponding distri- 
butions are simple Gaussians. 

However, sampling the spectral indices is a-priori not 
trivial — their distributions are highly non-Gaussian, and 
no analytical expressions exist. However, given that we 
already have run a MCMC analysis whose product is 
precisely the joint density F(scmb, ^d, /3s|d), the prob- 
lem is mostly solved. We could simply generate a full 
multi-dimensional histogram from the MCMC samples 
(for each pixel separately), and pick out the synchrotron 
index column that corresponds to the other currently 
fixed parameter values. Given this one-dimensional dis- 
tribution, we could then sample numerically using stan- 
dard techniques. 

For completeness, we note that adding more than 
one foreground component is a trivial extension of this 
scheme. Component amplitudes are added individually, 
while there is a choice for spectral indices — one may ei- 
ther sample these individually, as done above for the syn- 
chrotron index, or for greater efhciency, one may also ex- 
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WMAP 23 GHz band 




LFI 70 GHz band 




HFI 217 GHz band 




Fig. 2. — The "high-resolution" simulations used in this paper. 
Shown are the 23 GHz channel from the WMAP experiment, the 
70 GHz channel from the LFI experiment, and the 217 GHz channel 
from the HFI experiment. All maps are smoothed to a common 
resolution of 1° FWHM. 



ploit the multi-variate information given by the MCMC 
analysis. 

With the above prescription, it is finally possible to 
propagate the foreground uncertainties rigorously all the 
way from the observed data through to the CMB power 
spectrum, and therefore to the cosmological parame- 
ters. Further, with this approach one also obtains full- 
resolution sampled uncertainties of the component ampli- 
tudes, as opposed to the analytical approximations dis- 
cussed in the previous section, and a complete proba- 
bilistic description of the system is thereby established. 

In the form described above, the algorithm requires 
storage of a multi-dimensional histogram for each pixel. 
This is tremendously wasteful, and clearly not feasible 
due to memory limitations. A first necessary step is 
therefore to compress the likelihood information into a 
smaller set of numbers. One possible solution was pro- 



posed bv lVerde et al.l l|2003() . who fitted a 4th order poly- 
nomial to the log-likelihood, 

In £ = go + ^ + ^ q-i5i5j + 

t k v kl (^"^^ 

ijk ijkl 

Here qi are fit coefficients, and Si — (a^ — Q;^)/ai, be- 
ing parameter number i and its maximum likelihood 
value^^. With only 210 coefhcients for a six-parameter 
likelihood, this approximation at least satisfies the mem- 
ory requirements for practical usage in our problem. Fur- 
ther, it is straightforward to map out the required spec- 
tral index distributions given a set of component ampli- 
tudes using this form. Therefore, considering that the 
process parallelizes trivially, there seem to be no unsur- 
mountable computational problems connected to this ap- 
proach. 

While Gibbs sampling as currently implemented has 
problems with probing the low signal-to-noise ratio 
regime properly, it work s very well for signal-t o-noise ra- 
tios larger than unity (|Eriksen et al.l l2'004b|) . and this 
is exactly where the foreground uncertainties dominate. 
Therefore it seems reasonable to use the approach pre- 
sented here to analyze the high and intermediate signal- 
to-noise ratio regimes, propagating foreground uncertain- 
ties to th e final products, an d a standard "Master" -type 
analysis ijHivon et alJ 1200^ for the low signal-to-noise 
regime, at the cost of neglecting foreground uncertain- 
ties at these angular scales. 

6. EXAMPLE: APPLICATION TO PLANCK AND SIX- YEAR 
WMAP DATA 

We now apply the MCMC component separation 
method described in §0 to simulations of the current 
WMAP and the future Planck missions. We first give 
a detailed presentation of the simulations and data. We 
then study the behavior of the algorithm for one arbi- 
trarily chosen pixel, before considering the full sky map 
solutions. 

We point out that our main goal in this paper is to 
study the algorithm itself, and not to simulate an ac- 
tual data release. We therefore choose examples both 
with and without modeling errors, in order to illustrate 
problems that may be encountered in an analysis of real 
data. 

6.1. Simulations and models 

The simulations used in the following are constructed 
as a sum of a cosmological CMB signal, three fore- 
ground components (synchrotron, free-free, and ther- 
mal dust emission) and instrumental noise. We include 
five bands (centered at 23, 33, 41, 61 and 94 GHz) 
from WMAP, three bands (30, 44 and 70 GHz) from 
the Planck Low Frequency Instrument (LFI), and three 
bands (100, 143 and 217 GHz) from the High Frequency 
Instrument (HFI), for a total of eleven bands between 

Introducing normal parameters and/or higher-order expan- 
sions may be required for obtaining sufficient accuracy; see e.g., 
Chu, Kaplinghat & Knox (2003) for an application to cosmological 
parameters. 
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23 and 217 GHz^^. Specifications for each detector are 
given in Tabled 

CMB — The CMB component is assumed to be Gaus- 
sian distributed, with variances given by th e best- 
fit W MAP power-law power spectrum ( Benne tt et alJ 
I2003at iHmshaw et al.ll200l ISoergel eTallfeooaj) . includ- 
ing multipoles between £ = 2 and 1024. The signal 
reaHzation is filtered through the HEALPix pixel win- 
dow function and the instrument specific beam windows. 
(Since all the Planck beam windows are not available, we 
choose for simplicity to model even the WMAP beams 
as Gaussians with appropriate FWHM's.) 

Noise — The noise is assumed to be Gaussian and uncor- 
related, but non-uniform according to the scanning strat- 
egy of each detector. For WMAP, we assume a six- year 
mission, and rescale the published first-year sensitivity 
levels by 1/V6. For Planck, we adopt the requirement 
levels, which are a factor of two worse than the goals, for 
the baseline one-year mission. 

Synchrotron emission — The only all-sky map currently 
available to provide a template of Galactic synchrotron 
emission in HEALP ix format is the 408 MHz survey by 
iHaslam et al.l This has a resolution of only 51', 

thus additional power must be added on smaller angular 
scales for our purposes in this paper. We adopt the model 
of Giardino_et al. ( 2002), who estimated the amplitude 
and slope of the synchrotron angular power spectrum at 
low Galactic latitudes for I > 150. A Gaussian realization 
was then generated from that power spectrum to which 
was applied a Galactic modulation, multiplying t he sig- 
nal in each pixel by the ratio between the Hasla m et al.l 
()1982(1 template in that pixel and th e maximum. The fi- 
nal template is added to the original lHaslam et al.l l|1982fl 
map. An all-sky template for the synchrotron spectral 
index was estimated combining the all-sky data from 
iH aslam ct al. ( 1982), with the Northern sky observations 
of iReich & Reich] lll98(i[) at 1420 MHz and the southern 
sky counterpart of I Jonas et all 1)1998(1 at 2300 MHz". 
This constitutes our synchrotron model. 

Free-free emission — Current models of the free-free emis- 
sion exploit the expected correlation with Ha emission 
(see Dickinson et al. 2003 and references therein), and 
predict free-free emission (given in antenna temperature) 
as follows 



-,0.029- 
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Here y^Ha is the Ha amplitude, G = 
3MT°'^\v/vo^ff)-°-^'^ is the Gaunt factor ijFinkbeined 
I2OO3.') , Tc is the electron temperature, and i/q^s = 40 GHz 
is a reference frequency. 

A major uncertainty when using Ha as a template 
is due to the absorption of Ha by foreground dust, 

^® The three highest HFI frequency bands are not included in the 
analysis because they would introduce significant dust modeling 
errors; we simulate dust with a two-component model, but fit for 
a one-component model. 

^"^ This hybrid spectral index model should ultim ately be su- 
persed ed by the full-sky 1420 MHz survey described in lR.eich et alJ 



but this can be estim ated using the 100 /xm maps from 
iSchlegel et al.l 1(19981) . The largest uncertainty, however, 
is related to the fracti on of dust (f d) lying in front of the 
Ha-emitting region. iDickinson e t al. (2003) show that 
for regions 30° < I < 60° and 5° < \b\ < 15°, fd « 0.3, 
while for local high latitude regions such as Orion and 
the Gum Nebula, there is little or no absorption by dust 
{fd ^ 0). T he latter is su pported by the cross-correlation 
analysis of iBandav et a l. ( 2003) of the Ha data with 
C05i?-DMR, which contrasts with the value of ~ 0.5 
expected under the assumption that the ionized gas and 
dust are coextensive along the line of sight (i.e., uni- 
formly mixed) , as is assumed in the WMAP analysis of 
iBennett et all l(2003bD and iFinkbeiiie^ l(2004a() . 

We correct for dust absorption, by assuming a single 
component dust model, with a temperature of 18.3 K, 
and a absorption fraction of 0.33 up to a flux correspond- 
ing to 1 magnitude. We assume an electron temperature 
of Tc = 7000 K, and therefore an effective frequency scal- 
ing close to over the range of frequencies consid- 
ered here. For the future, a more accurate model, ac- 
counting for the steepening spectral index, could be im- 
plemented as a correction to the simple power-law model. 



Thermal dust — We adopt model 8 of iFinkbeiner et all 
(1999) for thermal dust emission, with parameters /i = 
0.0363, qi/q2 = 13, vo,d = 3000GHz, ai = 1.67, 
a2 = 2.70, Ti = 9.4 K, and T2 = 16.2 K (see equations 
1 ^11(1 . However, this is too many parameters to fit in- 
dividually, and we therefore adopt a simpler model for 
reconstruction (see equation EJ. Modeling errors of the 
sort to be expected with real data will result. 

Data processing — All simulations are initially made at a 
pixel resolution of iVgido = 512, corresponding to a pixel 
size of 7'. However, since our method requires identi- 
cal beam sizes for all frequency bands, we downgrade 
each band separately to 1° FWHM (determined by the 
52.8' FWHM beam of the 23 GHz WMAP band) and 
reduce the pixel resolution to A^sidc — 256 (by deconvolv- 
ing the original beam and pixel windows, and convolving 
the common 1° FWHM beam and lower resolution pixel 
window) . 

By downgrading the data, the noise specifications are 
also modified. To estimate the effective noise levels af- 
ter degradation, we therefore generate 1000 noise real- 
izations for each band, and downgrade these in the same 
manner as the actual data maps. The effective noise lev- 
els of the downgraded maps are then estimated by taking 
the standard deviation of the 1000 realizations. 

The data set described above constitutes our main 
simulation, and is referred to in the following as "high- 
resolution data". Examples are shown in Figure |21 For 
the non- linear parameter estimation step using MCMC, 
we noise levels must be lower, as discussed earlier. There- 
fore we smooth all maps with an additional 6° FWHM 
Gaussian beam, and downgrade the pixel resolution to 
A'^sido — 32. (This smoothing scale is not optimized, but 
it is sufficient for the purposes of the present paper.) 
Again, the effective noise levels are determined by Monte 
Carlo simulations. This smoothed data set is referred to 
as "low-resolution data" . 

Initial model map — Th e initial model map i s based on 
the WMAP KpO mask ((Bennett et alJl2003b(l . First the 
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Fig. 3. — Marginalized parameter probability distributions for an arbitrarily chosen pixel inside the Galactic plane, generated by MCMC 
as described in the text. The vertical lines show the true input value for the pixel. (The true value is not well defined for the dust spectral 
index, since dust is modeled by a two-component spectrum, while a one-component model is fitted. However, the stronger of the two dust 
components has an index parameter of oi = 1.67.) 



excluded region of the original mask is expanded by 10° 
in all directions. Then all accepted pixels (i.e., the high- 
latitude region) are assigned the model that includes 
CMB, synchrotron, and dust (both with free amplitude 
and spectral index), while the model for the rejected pix- 
els additionally includes a free- free amplitude. 

6.2. Results 

We now apply the method of §|Slto the simulated data 
set described above. 

6.2.1. Single pixel results 

We first examine the performance of the MCMC algo- 
rithm by studying one single pixel in the low-resolution 
data set, namely pixel number 6100, which is located in- 
side the Galactic plane at I = 58°, 6 = 0°. The reasons for 
choosing this pixel (or one like it) are twofold. First, the 
model for this pixel includes all three foreground compo- 
nents, and has thus a complicated probability structure. 
Second, the model is rejected by the goodness-of-fit test, 
and this example therefore illustrates the modeling error 
problem. 

As discussed earlier, the MCMC algorithm basically 
performs a random walk on the likelihood surface, pro- 
ducing a set of samples from which the likelihood may 
be estimated by constructing single or multi-dimensional 
histograms. Examples of such histograms are shown in 
Figures 01 and m 

The first figure shows the probability distributions for 
each of the six included parameters, marginalized over 
all other parameters. Comparing with the true input 
values (vertical lines), we see that the algorithm repro- 
duces the correct values, and also that the uncertainties 
are reasonable compared to the true errors. 



Second, in Figure 01 we show two-dimensional proba- 
bility distributions for the same parameters. The true 
values are marked by a box. Several points are worth 
noticing in this figure. First, all parameters are clearly 
correlated, and some specific pairs very tightly so. Exam- 
ples of the latter are dust amplitude versus dust spectral 
index, and synchrotron amplitude versus free-free ampli- 
tude. 

Second, many of the distributions are clearly non- 
Gaussian, and it is clear that a Gaussian approxima- 
tion at this stage will not yield reliable errors. Still, the 
structures appear to be reasonably well behaved, and in 
principle it may be possible to find analytical parame- 
ter transformations that could ease the computational 
burden. 

Third, while most of the true values lie inside the 3cr 
confidence regions, in one case, namely the synchrotron 
amplitude versus synchrotron index, it lies far outside the 
acceptable region. Another perspective of this is pro- 
vided by the x^j which for this pixel is 44. With five 
degrees of freedom (eleven frequencies and six free pa- 
rameters), this particular model is thus ruled out at the 
99.9999% confidence level. This is because we fit for a 
simpler model than the one used in the simulation: the 
data are smoothed by a wide 6° beam, and the thermal 
dust is fitted with a one-component model, whereas the 
simulation was based on a two-component model. This 
may also be seen Figure ^ where the fitted spectra for 
each component for this pixel are plotted. At low fre- 
quencies, the data points lie systematically above the 
fitted model, resulting in a clear rejection. 

However, even though the model is strongly rejected 
by the goodness-of-fit test, it is important to note all of 
the univariate distributions are still reasonable, and the 
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Fig. 4. — Marginalized two-dimensional probability distributions for the same pixel as in Figure|3] computed by MCMC. Boxes indicate 
the true input values, and the contours mark the peak and the 68, 95 and 99.7% confidence levels. 



CMB reconstruction is still useful. Therefore, a high 
does not necessarily imply that the pixel has to be dis- 
carded from further analysis, but rather that extra care 
has to be taken. Preferentially, the extra information in- 
dicated by the high ^ should be used to improve the 
model. 

6.2.2. Low-resolution full-sky maps 

We now consider the reconstructed full-sky maps, 
starting with the low-resolution maps as computed by 
the MCMC analysis. The individual component maps 
are shown in Figure |S1 with reconstructions given in the 
left column, differences between reconstructed and input 
maps in the middle column, and estimated errors in the 
right column. In the left column of Figure El we show 
the model map used in the analysis and the resulting 
goodness-of-fit distribution (as computed with param- 
eter means, not the maximum likelihood values). 

Starting with the goodness-of-fit map, we first note 
that we should expect ^ 13 at 2a confidence at high 
latitudes, since the model has six degrees of freedom in 
this region. This is indeed the case for two wide bands 
on each side of the Galactic plane, and both the model 
and the estimated parameters may therefore be accepted 



as they stand. However, at very high latitudes and, less 
surprisingly, at low latitudes, the goodness-of-fit is poor. 

In the left column of Figure [S] we show the six recon- 
structed parameter maps, and in the middle column the 
actual output versus input errors. Clearly, the method 
works very well, as the CMB sky map is virtually free of 
artifacts, with residuals less than 10% even in the inner 
Galactic plane. And with the exception of sharp bound- 
aries in the foreground reconstruction, due to different 
models used in different regions, the foreground results 
also look encouraging. 

However, as good as these results are, we warn the 
reader against interpreting them as an expected perfor- 
mance level for future missions. Even though our simula- 
tions are as realistic as possible given our current under- 
standing of foreground properties, they are certainly not 
as complicated as the real sky. Considerable modeling 
errors must be expected for real data sets, and sky cuts 
are very likely still required for future work. 

One note about the sharp boundaries seen in the fore- 
ground maps is in order. If the reconstructed maps are 
intended for foreground studies, such features are clearly 
not acceptable. In such cases, post-processing may be 



15 



Rcconslriiclioii Residual Estimated RMS error 




Fig. 5. — Low-resolution parameter maps reconstructed by MCMC. The left column shows the parameter estimates, the middle column 
shows the difference between reconstructed and input maps, and the right column shows the rms errors estimated by MCMC. From top to 
bottom we show: 1) the thermodynamic CMB temperature; 2) the synchrotron emission amplitude relative to 23 GHz; 3) the synchrotron 
spectral index; 4) the free-free emission amplitude relative to 33 GHz; 5) the thermal dust emission amplitude relative to 90 GHz; and 6) the 
thermal dust spectral index. 



required, for instance by smoothing the boundary by a 
Gaussian beam. On the other hand, if the maps are to be 
used for CMB power spectrum or cosmological parame- 
ter estimation, it is better to use the maps as they are, 



and propagate the pixel errors rehably; the boundaries 
are mainly due to different noise properties in the various 
regions. However, we point out that the distinct bound- 
aries seen in FigureElare at least partially due to a poorly 
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Fig. 6. — Example of iterative improvement of the model map. Two MCMC runs using the same data were made. The first run included 
CMB, synchrotron, and dust (both with free amplitude and spectral index) at high latitudes (green), and also free-free at low latitudes 
(red). The second run fixed the dust index at 1.55 at very high latitudes (blue). The bottom row shows the distribution in the two 
cases; notice the significant improvement at high latitudes resulting from removing a (non-critical) parameter from the system. Similar 
improvements at low latitudes could be made by trial-and-error. 



chosen model map, constructed from a WMAP galactic 
mask, rather than the specific simulation under consider- 
ation; manual tweaking would surely improve the results 
considerably. 

Returning for a moment to the goodness-of-fit map 
shown in Figure IHI and comparing with the rms maps 
shown in the right column of figure |S1 we see that the 
very high-latitude region with high corresponds di- 
rectly to the thermal dust spectral index map. Further, 
we also see that the dust amplitude is very low in the 
same region. The interpretation is clear: thermal dust 
is not well constrained in these regions because of its 
low amplitude, leading to poorly constrained spectral in- 
dices. This again propagates into the CMB component, 
and the overall fit is unacceptable. 

The solution to this problem seems obvious. Since the 
main problem is unconstrained dust spectral indices, we 
should manually fix them at some reasonable value. The 
potential bias introduced in the CMB and other compo- 
nents by doing so is very small because of the small dust 
amplitude found by the first analysis. We implement this 
by assigning a new model that fixes the dust spectral for 
all pixels with a dust spectral index rms larger than 0.15 
in the lower right panel of Figure |S1 The fixed spectral 
index value chosen is somewhat arbitrarily chosen to be 
1.55. The modified model map is shown in the top right 
panel of Figure El 

We now repeat the analysis, and obtain the goodness- 
of-fit map shown in the lower right panel of Figure IHI 
Clearly, introducing a new model at high latitudes had a 



very beneficial impact on the results. In principle, we 
could now proceed with similar considerations at low 
latitudes, and obtain reasonable fits over the full sky. 
However, since our main purpose in this paper is to illus- 
trate the method, we are content with the slightly revised 
model map shown in the top right panel of Figure |S1 and 
use this map in the rest of the paper. 

We now consider the error estimation accuracy of the 
MCMC algorithm. In the left panel of Figure [7| we plot 
a histogram of the relative CMB reconstruction error 
6 — (ATost — ATin)/crost, where ATest is the estimated 
CMB temperature, ATjn is the true value, and cJcst is 
the estimated error. If both the amplitude and the error 
are perfectly estimated, the pixel histogram will match a 
Gaussian distribution with vanishing mean and unit vari- 
ance. (In this plot, we include only pixels in the inter- 
mediate latitude region with a goodness-of-fit < 13.) 
Obviously, the algorithm works very well, as the bias is 
very small indeed and the estimated error is very close 
to the true error. 

Finally, in Figure |S1 we plot a histogram of the x^s 
of the same pixels, and compare it to a Xe distribution. 
Clearly, there is a small shift towards high values, indi- 
cating that there are pixels within the set for which the 
model is rejected. Once again, for an analysis of real 
data, we should then go back to the sky maps shown in 
Figures Inland El and try to locate the "offending" pixels. 

We conclude this section by making a few comments 
on the computational cost of the method. Running the 
MCMC analysis for each pixel is by far the most ex- 
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(a) Low-resolution MCMC results 



(b) High-resolution analytic results 



Fig. 7. — CMB reconstruction accuracy shown as a histogram of <5 = (ATeat — ATin)/(Test for pixels in the intermediate-latitude region. 
Results from the low-resolution MCMC analysis are shown in the left panel, and from the high-resolution analytical analysis in the right 
column. For perfect reconstruction of both CMB amplitude and error, both curves would match a Gaussian distribution with vanishing 
mean and unit variance (dashed curve). 
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Fig. 8. — Distribution of from the low-resolution MCMC anal- 
ysis, shown for the intermediate latitude region, with six degrees 
of freedom. The dashed curve shows the expected Xg distribution. 



pensive step of the algorithm. For weU-behaved pixels, 
we find that it takes on the order of 100 CPU seconds 
(divided over four processors per pixel) to reach the con- 
vergence criteria described above. For an Ngidc = 32 
map with 12,288 pixels, it therefore takes about 350 CPU 
hours per run. For clusters with of order 10^ processors, 
this is not a major problem. Further, since the algorithm 
scales with the number of pixels, and parallelizes trivially, 
it is not unreasonable to apply it at higher resolutions, 
say, at iVsido = 128 for 6000 CPU hours. 

6.2.3. High-resolution full-sky maps 

Having estimated the non-linear parameters by 
MCMC, the next step is to estimate the component am- 
plitudes from the full-resolution sky maps. As discussed 
earlier, this can be done cither with a Gibbs sampling ap- 
proach or with an analytic approach. In this paper, we 



choose the latter route, and leave the former to a future 
publication. 

The results from applying the method described in 
S I5.5I are shown in Figure El Once again, we see that 
the reconstructed parameter maps look visually com- 
pelling. There are few visible signs of contamination in 
the CMB reconstruction, and, indeed, even inside the 
central Galactic plane the errors are only a few tens of 

In the right panel of Figuredwe plot the relative CMB 
reconstruction error for the high-resolution map, as we 
did for the low-resolution map in the previous section. 
Two facts are clear from this plot. First, the bias is small, 
indicating that the analytic point estimator is quite accu- 
rate. Second, the histogram does not match the Gaussian 
distribution well, but is rather focused around smaller 
values. In other words, the errors are over-estimated by 
some small factor by the analytic error propagation for- 
mula. This should not be surprising, given the assump- 
tions that went into those calculations. Nevertheless, the 
estimated errors are in fact of the correct order of magni- 
tude, and they can therefore be used as a mental guide, 
although not for quantitative work. 

In Figure 1101 we plot the power spectrum of the re- 
constructed high-resolution CMB map, the true realiza- 
tion specific input spectrum, and the ensemble-averaged 
spectrum. The reconstructed spectrum was computed by 
full-sky integrals without noise weights. The results are 
therefore excellent — the reconstruction is virtually per- 
fect up to £ = 200, after which a small noise term starts 
to make an impact, before the 1° beam renders the recon- 
struction arbitrary at £ = 300. From this plot it seems 
clear that we were too conservative when choosing a 6° 
beam for the low-resolution analysis, and that the com- 
bined Planck and six-year WMAP data can easily handle 
higher resolutions. 




Fig. 9. — High-resolution parameter maps reconstructed by direct solution of linear systems, fixing spectral indices at low-resolution 
map values. The left column shows the parameter estimates, the middle column shows the difference between output and input maps, and 
the right column shows the analytically estimated errors. From top to bottom wc show: 1) the thermodynamic CMB temperature; 2) the 
synchrotron emission relative to 23 GHz; 3) the free-free emission relative to 33 GHz; and 4) the thermal dust emission relative to 90 GHz. 
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Fig. 10. — Reconstructed (blue curve) and input (red curve) 
CMB power spectra. The ensemble-averaged spectrum is shown as 
a black smooth curve. The reconstructed spectrum was computed 
by full-sky integration without noise weights or sky cut. 



7. OUTSTANDING PROBLEMS 

As demonstrated in the previous section, the compo- 
nent separation method presented in this paper works 
well on simulated data. However, there are a few out- 
standing issues that we have not taken into account, but 
that will have to be addressed prior to analysis of real 
data. The most important of these are gain and zero- 
point calibration, dipole corrections, noise correlations, 
and beam asymmetries. 

First, our algorithm requires all sky maps to be prop- 
erly calibrated with respect to gain and zero-point. Usu- 
ally, the gain is calibrated using the CMB dipole, but this 
can be difhcult for channels that are highly foreground- 
contaminated. Also, zero-point calibration is never easy. 

Closely related to these issues is dipole subtraction. 
The CMB dipole itself is hard to observe because of the 
large Doppler-dipole induced by the motion of the Solar 
system through space; it is usually subtracted in the map 
making process. Nevertheless, residual dipoles may cause 
serious problems for our algorithm unless accounted for. 

Finding reliable calibration methods for each of the 
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above problems is clearly essential. Fortunately, the 
number of degrees of freedom represented by these is- 
sues is quite small, and it may be possible to include 
them in the analysis by replacing the signal model 5*1, in 
equation El with 

1 £ 

S^ip) ^ g,,S^{p) ^ aim,,,Yim{p). (36) 

1=0 m=-l 

Following a rough calibration with external techniques, 
one could then use methods similar to those described in 
this paper to optimize the gains g^, and monopole and 
dipole coefficients a^rn.v 

From a conceptual point of view, correlated noise poses 
a more serious problem. For Planck, for example, the 
main effect will be to introduce stripes in the sky maps 
along the scanning path of the detectors, and locally, this 
has the same effect as an overall offset. Properly speak- 
ing, correlated noise is a problem for map making more 
than it is for component separation; however, residual 
effects can be expected. Only when actual data are in 
hand will it be clear how serious a problem it is. 

Finally, in this paper we have assumed that all detec- 
tors have identical beam response functions. This obvi- 
ously is not true for any real system, and corresponding 
errors are unavoidable. Fortunately, this is likely to have 
a negligible effect on the low-resolution analysis, since we 
smooth with an additional degree-scale beam, strongly 
suppressing small-scale asymmetries. Only in the high- 
resolution analysis is this effect likely to be important. 

8. DISCUSSION 

In this paper, we approach the problem of component 
separation with CMB data from the perspective of pa- 
rameter estimation. Our goal is to propagate foreground 
uncertainties all the way from observed data through to 
the final products, most importantly to the CMB power 
spectrum and cosmological parameters. This is more eas- 
ily facilitated with standard parameter techniques than 
with image processing techniques. 

We proposed and implemented one particular algo- 
rithm for performing this task, based on multi-frequency 
parametric model fits established by means of a hybrid of 
Markov Chain Monte Carlo and analytic methods. The 
method was then shown to work very well on simulated 
data, with properties corresponding to those of the future 
Planck and six-year WMAP experiments. 

We also outlined how to propagate the foreground- 
induced errors to the CMB power spectrum and cos- 

^® http:/ /www. eso.org/science/healpix/ 



mological parameter errors, using the output from the 
MCMC analysis presented here as the input in a Gibbs 
sampling algorithm. As always, only an actual imple- 
mentation will prove whether this method works or not, 
but the theoretical groundwork appears to be sound, and 
no unsurmountable computational problems have been 
identified. Therefore, if this approach proves success- 
ful, we will have a complete, mathematically consistent, 
end-to-end solution to the foreground problem in CMB 
analysis. 

While we only considered temperature anisotropy ob- 
servations in the present paper, the method is completely 
general, and can equally well handle polarization mea- 
surements, as will be demonstrated in a future study. 
We will also apply the method to two specific problems. 
First, we will study the optimization of frequency cov- 
erage and signal-to-noise ratio in future polarization ex- 
periments. Since our method provides error bars on all 
estimated quantities, it is straightforward to compare dif- 
ferent experiment designs. Modeling errors will be an 
integral part of this work, since such uncertainties have 
a direct impact on the optimal frequency range to ob- 
serve. Second, we will apply the method to the currently 
available WMAP data. 
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